Astronomical Image Calibration

Tom Sherlock

User Interfaces, Wolfram Research, Inc.

What is Image Calibration?

Details in astronomical objects, apart from stars and planets, are generally quite faint and hard to see, even in long-exposure images.

If you do scientific imaging, like photometry or spectroscopy, handling image noise is quite important.

In order to maximize signal and minimize noise in astronomical images, several sources of noise must be removed, through image calibration.  Mathematica’s statistical and image processing capabilities make it ideal for this purpose.

Signal to Noise Ratio (SNR)

Combining many exposures maximizes the SNR

SNR is proportional to the square root of the number of frames

The longer each frame’s exposure, the better because we will have more photons from the object in each CCD pixel and the signal from those photons can be more easily distinguised from the noise coming from a variety of sources.

Major Sources of Noise

Thermal effects on the CCD chip (partial solution: cool the chip)

Read noise (uncertainty in the number of actual electrons read out of a CCD pixel)

Optical effects (dirt on the chip, or uneven illumination, for example).

Skyglow (solution: go to a darker site, or use narrow-band filters)

Minimizing these effects involves collecting additional calibration data besides the actual photons collected from the stars.

Different Types of Calibration Data

Dark Frames

These are exposures used to approximate the noise associated with thermal effects of the CCD camera’s circuitry.  They are taken with the camera covered by a shutter (or a hat) and don’t include photons from the sky. These will also include ‘hot’ pixels (defects in the CCD array) which need to be subtracted out. 

For this reason, it is important that the CCD has cooled down to its ‘working’ temperature when these are taken.  Better cameras allow you to set the temperature and thus keep a library of dark frames.

It is also important that you take about the same  number of dark frames as normal light frames and have their exposures match, so that you may closely approximate the amount of noise in the actual images.

Bias Frames

These are extremely short, dark, exposures used to gauge the noise associated with reading data from the CCD chip’s pixels.

Flat Frames

These are short, evenly illuminated light exposures used to remove the effect of uneven illumination on the CCD chip, caused, for example, by light falloff or dust particles.  

Flat frames become more important with larger CCD chips which may exceed the 100% illuminated field of the telescope.

Flat frames are generally exposed by aiming the scope at an evenly illuminated surface, like a special calibration screen, or even the twilit sky.

This is an approximation of the effect we try to eliminate by taking a flat frame:

ImageMultiply[,]

Light Frames

These are the actual images taken of the astronomical object of interest.

Calibration Pipeline

The calibration process consists of building a set of “calibration frames” and then applying them to the collected light frames before stacking.

Stacking Pipeline

Stacking combines aligned frames through the use of the Mean and Median functions in Mathematica to create a composite frame with a larger signal-to-noise ratio.

Example: M34 (Spiral Nebula)

20 x 60 second guided subframes.

Image Stacking Utilities

Create a function to take a list of images and create a new image by taking the Median or Mean of all
the list elements.

Stack using Median

ImageMedianList[sepList_List]:=
Module[{},
ColorCombine[{Image[Median[Map[ImageData,Map[#[[1]]&,sepList]]]],
Image[Median[Map[ImageData,Map[#[[2]]&,sepList]]]],
Image[Median[Map[ImageData,Map[#[[3]]&,sepList]]]]}]]

ImageMedianFlat[imageList_List]:=
Module[{sepList},
sepList=Map[ColorSeparate,imageList];
ColorCombine[{Image[Median[Map[ImageData,Map[#[[1]]&,sepList]]]],
Image[Median[Map[ImageData,Map[#[[2]]&,sepList]]]],
Image[Median[Map[ImageData,Map[#[[3]]&,sepList]]]]}]]

ImageMedian[imageList_List]:=If[Head[imageList[[1]]]==List,ImageMedianList[imageList],ImageMedianFlat[imageList]]

Stack using Mean

ImageMeanList[sepList_List]:=
Module[{},
ColorCombine[{Image[Mean[Map[ImageData,Map[#[[1]]&,sepList]]]],
Image[Mean[Map[ImageData,Map[#[[2]]&,sepList]]]],
Image[Mean[Map[ImageData,Map[#[[3]]&,sepList]]]]}]]

ImageMeanFlat[imageList_List]:=
Module[{sepList},
sepList=Map[ColorSeparate,imageList];
ColorCombine[{Image[Mean[Map[ImageData,Map[#[[1]]&,sepList]]]],
Image[Mean[Map[ImageData,Map[#[[2]]&,sepList]]]],
Image[Mean[Map[ImageData,Map[#[[3]]&,sepList]]]]}]]

ImageMean[imageList_List]:=If[Head[imageList[[1]]]==List,ImageMeanList[imageList],ImageMeanFlat[imageList]]

Stack using Mean after dropping the fraction f of the smallest and largest elements

TrimmedImageMeanList[sepList_List, f_]:=
Module[{},
ColorCombine[{Image[TrimmedMean[Map[ImageData,Map[#[[1]]&,sepList]],f]],
Image[TrimmedMean[Map[ImageData,Map[#[[2]]&,sepList]],f]],
Image[TrimmedMean[Map[ImageData,Map[#[[3]]&,sepList]],f]]}]]

TrimmedImageMeanFlat[imageList_List, f_]:=
Module[{sepList},
sepList=Map[ColorSeparate,imageList];
ColorCombine[{Image[TrimmedMean[Map[ImageData,Map[#[[1]]&,sepList]],f]],
Image[TrimmedMean[Map[ImageData,Map[#[[2]]&,sepList]],f]],
Image[TrimmedMean[Map[ImageData,Map[#[[3]]&,sepList]],f]]}]]

TrimmedImageMean[imageList_List, f_]:=If[Head[imageList[[1]]]==List,TrimmedImageMeanList[imageList,f],TrimmedImageMeanFlat[imageList,f]]

CreateMasterBiasFrame[biasFrames_List]:=ImageMedian[biasFrames];

CreateMasterDarkFrame[darkFrames_List, masterBiasFrame_]:=
ImageSubtract[ImageMedian[darkFrames],masterBiasFrame];

CreateMasterFlatFrame[flatFrames_List, masterDarkFrame_]:=
ImageSubtract[ImageMedian[flatFrames],masterDarkFrame];

CreateMasterFrames[biasFrames_List, darkFrames_List, flatFrames_List]:=
Module[{mBF,mDF,mFF},
mBF=CreateMasterBiasFrame[biasFrames];
mDF=CreateMasterDarkFrame[darkFrames, mBF];
mFF=CreateMasterFlatFrame[flatFrames, mDF];
List[mBF, mDF, mFF]
]

CreateMasterFrames[biasFrames_List, darkFrames_List]:=
Module[{mBF,mDF},
mBF=CreateMasterBiasFrame[biasFrames];
mDF=CreateMasterDarkFrame[darkFrames, mBF];
List[mBF, mDF]
]

CreateMasterFrames[biasFrames_List]:=
Module[{mBF},
mBF=CreateMasterBiasFrame[biasFrames];
List[mBF]
]

CalibrateImage[lightFrame_, masterBiasFrame_, masterDarkFrame_, masterFlatFrame_]:=ImageMultiply[ImageSubtract[ImageSubtract[lightFrame,masterBiasFrame],masterDarkFrame],ColorNegate[masterFlatFrame]];

CalibrateImage[lightFrame_, masterBiasFrame_, masterDarkFrame_]:=ImageSubtract[ImageSubtract[lightFrame,masterBiasFrame],masterDarkFrame];

CalibrateImage[lightFrame_, masterDarkFrame_]:=ImageSubtract[lightFrame,masterDarkFrame];

Processing Calibration Data

Bring in the FITS files for the different frames.

dataLocation="c:\\users\\sherlock\\desktop\\2014WTC\\ImageData\\";
biasColorFrames=Map[Import,FileNames[dataLocation<>"2014_10_8_M34_Bias_*.fit"]];
darkColorFrames=Map[Import,FileNames[dataLocation<>"2014_10_8_M34_Dark*.fit"]];
lightColorFrames=Map[Import,FileNames[dataLocation<>"2014_10_8_M34_Light2_*.fit"]];

Build the ‘master’ dark and bias frames by taking the median of all the dark and bias frames.  In this example, I don’t have any flat frames, so they aren’t used.

Here CreateMasterFrames takes the Mean of the sets of bias and dark frames to create a list of calibration frames.

masterFrames=CreateMasterFrames[biasColorFrames, darkColorFrames];

Perform the calibration pipeline on our image data to subtract out the noise associated with thermal and read effects.

calibratedFrames=Map[CalibrateImage[#, masterFrames[[1]], masterFrames[[2]]]&,Map[ColorCombine, lightColorFrames]];

Aligning Frames

We need to align the calibrated frames in preparation for stacking.

Stacking reduces the random noise by statistically averaging it away.

Mathematica can automatically align the frames, but we can generally get better results by manually aligning them.

Fortunately, it is easy to use Mathematica to produce a user interface to align the frames.

Frame Alignment Process

Locate a non-saturated star image in each frame.
The coordinates of this star will be used to compute the transforms that can
be used to shift the frames onto one another

myPoints=Table[{1,1},{Length[calibratedFrames]}];

Manipulate[LocatorPane[Dynamic[myPoints[[i]]],calibratedFrames[[i]]],{i,1,Length[calibratedFrames],1}]

We don’t have to use the first frame as the key frame, we can pick any frame, but
here we will just use the first one.

keyFrame=myPoints[[1]];
geoTransforms=FindGeometricTransform[{keyFrame},{#}]&/@myPoints;
transformImage[im_,i_]:=ImagePerspectiveTransformation[im,geoTransforms[[i]][[2]],DataRange->Full]
transIm={};

Build the list of transformed image frames

For[i=1,i<Length[geoTransforms],i++,AppendTo[transIm,transformImage[calibratedFrames[[i]],i]]];

Stacking Frames

Now that we have them all aligned, we can use our statistical functions, Median and Mean, to obtain averaged images, one using Mean and one using Median.

stackedImMedian=Image[Median[Map[ImageData,transIm]]]

stackedImageMean=Image[Mean[Map[ImageData,transIm]]]

Final Image Processing

m34Adj=ImageAdjust[stackedImMedian,{.4,1}]

Sharpen[m34Adj]

Conclusion

Mathematica’s powerful language, statistical, image processing and user interface features make it the ideal tool for performing the complex tasks needed to reduce astronomical data.